library(spatstat)
library(sf)
library(sf)

Data Exploration

Load in our Covariate Data:

load('BC_Covariates.Rda')
BC_Cov <- DATA
BC_Cov <- na.omit(BC_Cov)
ls()
## [1] "BC_Cov" "DATA"
head(summary(BC_Cov))
##            Length Class           Mode
## Window      1     SpatialPolygons S4  
## Elevation  10     im              list
## Forest     10     im              list
## HFI        10     im              list
## Dist_Water 10     im              list

Load in our species data and convert it to a ppp object:

bees <- read.csv('0014081-250402121839773.csv', sep = "\t")
coordinates(bees) <- ~decimalLongitude + decimalLatitude
proj4string(bees) <- CRS("+proj=longlat +datum=WGS84")
aea_crs <- CRS("+proj=aea +lat_0=45 +lon_0=-126 +lat_1=50 +lat_2=58.5 +x_0=1000000 +y_0=0 +datum=NAD83 +units=m +no_defs")
bees_proj <- spTransform(bees, aea_crs)
bee_coords <- coordinates(bees_proj)
head(bees,2)
##         coordinates    gbifID                           datasetKey occurrenceID
## 1 (-118.9161, 50.4) 735215733 275319e1-f91c-406f-b239-62cb9d4185cb   LEM0258825
## 2 (-118.9161, 50.4) 735215732 275319e1-f91c-406f-b239-62cb9d4185cb   LEM0258824
##    kingdom     phylum   class       order   family   genus              species
## 1 Animalia Arthropoda Insecta Hymenoptera Vespidae Vespula Vespula pensylvanica
## 2 Animalia Arthropoda Insecta Hymenoptera Vespidae Vespula Vespula pensylvanica
##   infraspecificEpithet taxonRank                           scientificName
## 1                   NA   SPECIES Vespula pensylvanica (de Saussure, 1857)
## 2                   NA   SPECIES Vespula pensylvanica (de Saussure, 1857)
##   verbatimScientificName verbatimScientificNameAuthorship countryCode
## 1                                        (Saussure, 1857)          CA
## 2                                        (Saussure, 1857)          CA
##         locality    stateProvince occurrenceStatus individualCount
## 1 Trinity Valley British Columbia          PRESENT               1
## 2 Trinity Valley British Columbia          PRESENT               1
##                       publishingOrgKey coordinateUncertaintyInMeters
## 1 182287da-fe0f-457f-ac06-2b8f4c7680fd                            NA
## 2 182287da-fe0f-457f-ac06-2b8f4c7680fd                            NA
##   coordinatePrecision elevation elevationAccuracy depth depthAccuracy
## 1                  NA        NA                NA    NA            NA
## 2                  NA        NA                NA    NA            NA
##    eventDate day month year taxonKey speciesKey      basisOfRecord
## 1 1939-04-28  28     4 1939  1311698    1311698 PRESERVED_SPECIMEN
## 2 1939-05-02   2     5 1939  1311698    1311698 PRESERVED_SPECIMEN
##     institutionCode collectionCode catalogNumber recordNumber identifiedBy
## 1 McGill University           LEMQ    LEM0258825           NA             
## 2 McGill University           LEMQ    LEM0258824           NA             
##   dateIdentified license      rightsHolder  recordedBy typeStatus
## 1                CC0_1_0 McGill University C.V. Morgan         NA
## 2                CC0_1_0 McGill University C.V. Morgan         NA
##   establishmentMeans          lastInterpreted mediaType
## 1                 NA 2025-02-01T03:52:23.555Z          
## 2                 NA 2025-02-01T03:52:23.565Z          
##                                                                                                  issue
## 1 GEODETIC_DATUM_ASSUMED_WGS84;OCCURRENCE_STATUS_INFERRED_FROM_INDIVIDUAL_COUNT;INSTITUTION_MATCH_NONE
## 2 GEODETIC_DATUM_ASSUMED_WGS84;OCCURRENCE_STATUS_INFERRED_FROM_INDIVIDUAL_COUNT;INSTITUTION_MATCH_NONE
window_sf <- st_as_sf(BC_Cov$Window)
win <- as.owin(window_sf) 

bee_ppp <- ppp(x = bee_coords[,1], y = bee_coords[,2], window = win)
table(inside.owin(bee_ppp$x, bee_ppp$y, win))  
## 
## TRUE 
## 1112
#summary(bee_ppp)

We can plot the elevation as a persp figure, with our wasp data points overlayed on top:

library(viridis)

fig <- persp(BC_Cov$Elevation,
             theta = -15, phi = 30, # rotation
             expand = 7, # z-axis expansion
             border = NA, #remove grid borders
             #apron = TRUE, #apron around edge
             shade = 0.3, # shading
             box = FALSE, # axes on/off
             main = "", # title
             visible = TRUE, #Supporting calculations
             colmap = viridis(200)) # colour pallet

perspPoints(bee_ppp,
            Z = BC_Cov$Elevation, 
            M = fig,
            pch = 20, 
            cex = 0.6)

Additionally, we can plot each covariate with the points overlayed on top:

plot(BC_Cov$Elevation, main = "Elevation + Bee Observations")
plot(bee_ppp, add = TRUE, pch = 20, cex = 0.6)

plot(BC_Cov$Forest, main = "Forest Cover + Bee Observations")
plot(bee_ppp, add = TRUE, pch = 20, cex = 0.6)

plot(BC_Cov$HFI, main = "Human Footprint Index + Bee Observations")
plot(bee_ppp, add = TRUE, pch = 20, cex = 0.6)

plot(BC_Cov$Dist_Water, main = "Distance to Water + Bee Observations")
plot(bee_ppp, add = TRUE, pch = 20, cex = 0.6)

We can also plot the density of each covariate at each bee obserevation compared to the density of each covariate across the sampling window. A difference in the curves may tell us that Yellowjackets may prefer or avoid areas with certain values of that covariate.

library(ggplot2)

#elevation 
elev_df <- data.frame(
  value = c(BC_Cov$Elevation[], BC_Cov$Elevation[bee_ppp]),
  group = c(rep("Background", length(BC_Cov$Elevation[])),
            rep("Bee Points", length(BC_Cov$Elevation[bee_ppp])))
)

ggplot(elev_df, aes(x = value, fill = group)) +
  geom_density(alpha = 0.4) +
  labs(title = "Elevation Density", x = "Elevation", y = "Density") +
  scale_fill_manual(values = c("blue", "red")) +
  theme_minimal()

#forest cover
forest_df <- data.frame(
  value = c(BC_Cov$Forest[], BC_Cov$Forest[bee_ppp]),
  group = c(rep("Background", length(BC_Cov$Forest[])),
            rep("Bee Points", length(BC_Cov$Forest[bee_ppp])))
)

ggplot(forest_df, aes(x = value, fill = group)) +
  geom_density(alpha = 0.4) +
  labs(title = "Forest Cover Density", x = "Forest", y = "Density") +
  scale_fill_manual(values = c("blue", "red")) +
  theme_minimal()

#distance to water
water_df <- data.frame(
  value = c(BC_Cov$Dist_Water[], BC_Cov$Dist_Water[bee_ppp]),
  group = c(rep("Background", length(BC_Cov$Dist_Water[])),
            rep("Bee Points", length(BC_Cov$Dist_Water[bee_ppp])))
)

ggplot(water_df, aes(x = value, fill = group)) +
  geom_density(alpha = 0.4) +
  labs(title = "Distance to Water Density", x = "Dist_Water", y = "Density") +
  scale_fill_manual(values = c("blue", "red")) +
  theme_minimal()

hfi_df <- data.frame(
  value = c(BC_Cov$HFI[], BC_Cov$HFI[bee_ppp]),
  group = c(rep("Background", length(BC_Cov$HFI[])),
            rep("Bee Points", length(BC_Cov$HFI[bee_ppp])))
)

ggplot(hfi_df, aes(x = value, fill = group)) +
  geom_density(alpha = 0.4) +
  labs(title = "HFI Density", x = "HFI", y = "Density") +
  scale_fill_manual(values = c("blue", "red")) +
  theme_minimal()

First Moment Descriptive Statistics

bee_win_km <- rescale(Window(bee_ppp), s = 1000, unitname = "km")
npoints(bee_ppp) / area(bee_win_km)
## [1] 0.001172671

In this section, we estimated the first-order intensity of bee observations across the study region. Initially, the intensity was calculated in units of points per square meter, resulting in an extremely small value (~1.17e-09), which is difficult to interpret. To make the results more intuitive, we rescaled the spatial window to units of square kilometers, yielding a new intensity estimate of ~0.00117 bees/km².

This very low value suggests that bee observations are extremely sparse over the study area. While this could indicate a naturally low density of bees in the region, it is also possible that the bees are concentrated in specific locations, rather than being evenly distributed. This motivates further analysis of spatial inhomogeneity, using methods like quadrat counts, kernel density estimation, and second-order statistics such as Ripley’s K-function, which will be explored in the subsequent sections.

library(spatstat.explore)

Q <- quadratcount(bee_ppp, nx = 5, ny = 5)
plot(bee_ppp, pch = 20, cex = 0.5)
## Warning in plot.ppp(bee_ppp, pch = 20, cex = 0.5): 47 illegal points also
## plotted
plot(Q, add = TRUE)

quadrat.test(bee_ppp, nx = 5, ny = 5)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  bee_ppp
## X2 = 7290.7, df = 20, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 21 tiles (irregular windows)
plot(intensity(Q, image = TRUE),
     main = "Vespula pensylvanica intensity")

plot(bee_ppp,
     pch = 16,
     cex = 0.6,
     cols = "white",
     add = TRUE)
## Warning in plot.ppp(bee_ppp, pch = 16, cex = 0.6, cols = "white", add = TRUE):
## 47 illegal points also plotted
plot(bee_ppp,
     pch = 16,
     cex = 0.5,
     cols = "black",
     add = TRUE)
## Warning in plot.ppp(bee_ppp, pch = 16, cex = 0.5, cols = "black", add = TRUE):
## 47 illegal points also plotted

The quadrat-based intensity map of Vespula pensylvanica shows clear spatial heterogeneity in the species’ distribution. The highest intensities are concentrated in the southern and southeastern parts of the study region, while most of the northern areas show near-zero intensity. This suggests that V. pensylvanica is either ecologically constrained to certain environments or is under-sampled in other areas.

Interpretation

In this analysis, I performed a quadrat test to evaluate whether the spatial distribution of wasp observations in British Columbia follows a pattern of complete spatial randomness (CSR). The region was divided into a 5x5 grid (with 21 valid quadrats), and bee counts were tallied within each quadrat.

The resulting map shows clear spatial heterogeneity in wasp observations. Some quadrats, especially in the southern region, have very high counts (e.g., 322 and 66), while others—particularly in the north—have counts close to or exactly zero. This visual pattern already suggests a non-random distribution.

The quadrat test produced the following results:
Chi-squared = 7290.7, df = 20, p-value < 2.2e-16
These results indicate a highly significant deviation from CSR. The extremely low p-value allows us to reject the null hypothesis that bee occurrences are randomly distributed across space.

Conclusion

There is strong evidence that the bees are not distributed randomly, but rather show significant spatial clustering. This may reflect environmental, ecological, or human factors influencing wasp presence in the region.

library(spatstat.geom)
library(spatstat.explore)

lambda_u_hat <- density(bee_ppp) 

plot(lambda_u_hat,
     main = "Kernel estimate of bee_ppp intensity")

plot(bee_ppp,
     pch = 16,
     cex = 0.6,
     cols = "white",
     add = TRUE)
## Warning in plot.ppp(bee_ppp, pch = 16, cex = 0.6, cols = "white", add = TRUE):
## 47 illegal points also plotted
plot(bee_ppp,
     pch = 16,
     cex = 0.5,
     cols = "black",
     add = TRUE)
## Warning in plot.ppp(bee_ppp, pch = 16, cex = 0.5, cols = "black", add = TRUE):
## 47 illegal points also plotted

The kernel density estimate of bee_ppp shows that wasp observations are strongly clustered in the southern part of the study area. The intensity surface highlights areas with the highest densities in yellow, suggesting these are potential hotspots for wasp activity or observation effort. The presence of many zero- or low-intensity zones in the north may indicate ecological constraints, lower sampling effort, or habitat unsuitability.

library(spatstat.geom)
library(spatstat.explore)

par(mfrow = c(1, 3), mar = rep(0.1, 4)) 

# 1. Diggle bandwidth
plot(density(bee_ppp, sigma = bw.diggle(bee_ppp)), 
     ribbon = FALSE,
     main = "Bandwidth: Diggle")

# 2. PPL bandwidth
plot(density(bee_ppp, sigma = bw.ppl(bee_ppp)), 
     ribbon = FALSE,
     main = "Bandwidth: PPL")

# 3. Adaptive KDE
lambda_u_hat_adaptive <- adaptive.density(bee_ppp, method = "kernel")
plot(lambda_u_hat_adaptive,
     main = "Adaptive kernel estimate")

The three kernel density maps illustrate different smoothing strategies for the bee point pattern (bee_ppp). The Diggle and PPL methods use fixed bandwidths, with PPL producing slightly sharper intensity contrasts. The adaptive kernel method adjusts the bandwidth based on local point density, resulting in enhanced resolution in dense areas while maintaining smoother estimates in sparse regions. However, all three maps show similar general patterns, with high-intensity hotspots concentrated in the southern region of the study area.

library(spatstat.explore)

R <- bw.ppl(bee_ppp)

LR <- scanLRTS(bee_ppp, r = R)

plot(LR, main = "Local likelihood ratio test (bee_ppp)")

contour(LR, add = TRUE)
plot(bee_ppp, add = TRUE, pch = 16, cex = 0.4, cols = "black")
## Warning in plot.ppp(bee_ppp, add = TRUE, pch = 16, cex = 0.4, cols = "black"):
## 47 illegal points also plotted

pvals <- eval.im(pchisq(LR,
                        df = 1,
                        lower.tail = FALSE))
plot(pvals, main = "Local p-values")

The local likelihood ratio test (LRT) highlights strong spatial clustering of wasp observations in the southern region. The corresponding p-value map confirms that these clusters are statistically significant, with p-values close to zero (blue areas). This suggests that the observed clustering is unlikely to have occurred under complete spatial randomness (CSR) and indicates potential environmental or ecological drivers influencing bee distribution.

bee_elev <- BC_Cov$Elevation
elev <- bee_elev

b <- quantile(elev, probs = (0:4)/4, type = 2)
Zcut <- cut(elev, breaks = b)
V <- tess(image = Zcut)

quadratcount(bee_ppp, tess = V)
## Warning in quadratcount.ppp(bee_ppp, tess = V): Tessellation does not contain
## all the points of X
## tile
##          (-130,761]       (761,1.1e+03]  (1.1e+03,1.46e+03] (1.46e+03,3.56e+03] 
##                1040                  54                   9                   8
rho <- rhohat(bee_ppp, elev)
## Warning: Values for 1 query point lying outside the pixel image domain were
## estimated by projection to the nearest pixel
plot(rho, main = "Intensity of wasps vs elevation")

library(spatstat.explore)

env <- envelope(bee_ppp, Gest, nsim = 39)
## Generating 39 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
## 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 
## 39.
## 
## Done.
plot(env, main = "CSR test via G-function envelope")

A spatial point pattern analysis was conducted to assess the relationship between bee intensity and elevation, as well as deviations from complete spatial randomness (CSR). The non-parametric estimate showed that bee intensity is strongly concentrated at lower elevations. A Monte Carlo envelope test using the G-function further revealed that the observed bee locations exhibit significant clustering, with the observed G-function lying well above the simulated CSR envelopes.

Second Moment Descriptive Statistics

Morisita’s Index plot

rect_win <- as.owin(as.rectangle(win))
bee_ppp_rect <- ppp(x = bee_coords[,1], y = bee_coords[,2], window = rect_win)
## Warning: data contain duplicated points
miplot(bee_ppp_rect, main = "Morisita Index (Rectangular Window)", pch = 16, col = "#046C9A")

Ripley’s K-function estimate

xrange <- range(bee_coords[,1])
yrange <- range(bee_coords[,2])
win_rect <- owin(xrange = xrange, yrange = yrange)

unique_coords <- bee_coords[!duplicated(bee_coords), ]
bee_ppp_unique <- ppp(x = unique_coords[,1], y = unique_coords[,2], window = win_rect)
k_bee <- Kest(bee_ppp_unique)
plot(k_bee, main = "Ripley’s K-function (No Duplicates)", lwd = 2)

K-function with bootstrapped confidence intervals (α=0.05)

E_bee <- envelope(bee_ppp, Kest, correction = "border", rank = 1, nsim = 19, fix.n = TRUE)
## Generating 19 simulations of CSR with fixed number of points  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 
## 19.
## 
## Done.
plot(E_bee, main = "", lwd = 2)

K-function with bootstrapped confidence intervals (α=0.01)

E_bee_99 <- envelope(bee_ppp, Kest, correction = "border", rank = 1, nsim = 99, fix.n = TRUE)
## Generating 99 simulations of CSR with fixed number of points  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
## 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
## 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
## 99.
## 
## Done.
plot(E_bee_99, main = "", lwd = 2)

Estimate intensity for inhomogeneous K and PCF

lambda_bee <- density(bee_ppp, sigma = bw.ppl(bee_ppp), positive = TRUE)
cat("Mean estimated intensity (λ):", mean(lambda_bee$v, na.rm = TRUE), "\n")
## Mean estimated intensity (λ): 9.824262e-10
cat("Max estimated intensity (λ):", max(lambda_bee$v, na.rm = TRUE), "\n")
## Max estimated intensity (λ): 1.962303e-07

Inhomogeneous K-function

Kinhom_bee <- Kinhom(bee_ppp, lambda = lambda_bee)
## number of data points exceeds 1000 - computing border correction estimate only
plot(Kinhom_bee, theo ~ r, main = "", col = "grey70", lty = "dashed", lwd = 2)
## Warning: Internal error: fvlabels truncated the function name
## Warning: Internal error: fvlabels truncated the function name
## Warning: Internal error: fvlabels truncated the function name
plot(Kinhom_bee, border ~ r, col = "#046C9A", lwd = 2, add = T)
## Warning: Internal error: fvlabels truncated the function name
## Warning: Internal error: fvlabels truncated the function name
## Warning: Internal error: fvlabels truncated the function name

Simulation envelope for inhomogeneous K-function

E_bee_inhom <- envelope(bee_ppp, Kinhom, simulate = expression(rpoispp(lambda_bee)), correction = "border", rank = 1, nsim = 19, fix.n = TRUE)
## Warning in envelope.ppp(bee_ppp, Kinhom, simulate =
## expression(rpoispp(lambda_bee)), : fix.n and fix.marks were ignored, because
## 'simulate' was given
## Generating 19 simulations by evaluating expression  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 
## 19.
## 
## Done.
par(mfrow = c(1, 2))
plot(E_bee_inhom, main = "", lwd = 2)
plot(E_bee_inhom, xlim = c(0, 50000), main = "", lwd = 2)

Fitting and Validating Poisson Point Process Models

Finding out the shape of each Covariates in relation to the density of the species in BC using Relative Density Plots.

elev <- BC_Cov$Elevation
forest <- BC_Cov$Forest
hfi <- BC_Cov$HFI
distwater <- BC_Cov$Dist_Water

elev[is.na(elev)] <- median(elev, na.rm = TRUE)
forest[is.na(forest)] <- median(forest, na.rm = TRUE)
hfi[is.na(hfi)] <- median(hfi, na.rm = TRUE)
distwater[is.na(distwater)] <- median(distwater, na.rm = TRUE)

rho_elev <- rhohat(bee_ppp, elev)
## Warning: Values for 1 query point lying outside the pixel image domain were
## estimated by projection to the nearest pixel
rho_forest <- rhohat(bee_ppp, forest)
## Warning: Values for 1 query point lying outside the pixel image domain were
## estimated by projection to the nearest pixel
rho_hfi <- rhohat(bee_ppp, hfi)
## Warning: Values for 1 query point lying outside the pixel image domain were
## estimated by projection to the nearest pixel
rho_distwater <- rhohat(bee_ppp, distwater)
plot(rho_elev,
     main = "",
     xlab = "Elevation",
     xlim = c(0, max(elev)))

plot(rho_forest,
     main = "",
     xlab = "Forest Cover"
     )

plot(rho_hfi,
     main = "",
     xlab = "HFI"
     )

plot(rho_distwater,
     main = "",
     xlab = "Distance Water"
     )

# Check for collinearity
cor.im(elev, forest, hfi, distwater, use = "pairwise.complete.obs")
##             ..1         ..2         ..3         ..4
## ..1  1.00000000 -0.26204718 -0.26625626 -0.03426387
## ..2 -0.26204718  1.00000000  0.06615541  0.04825439
## ..3 -0.26625626  0.06615541  1.00000000  0.13229408
## ..4 -0.03426387  0.04825439  0.13229408  1.00000000

From the plots, we propose the following model with a combination of linear and quadratic terms

\[ \log(\lambda) = \beta_0 + \beta_1 \cdot \text{elev} + \beta_2 \cdot \text{forest} + \beta_3 \cdot \text{forest}^2 + \beta_4 \cdot \text{hfi} + \beta_5 \cdot \text{hfi}^2 + \beta_6 \cdot \text{distwater} \]

We then need to evaluate the quality of our model (how well the model is able to capture the variability in the data)

covariates_list <- list(elev = elev, forest = forest, hfi = hfi, distwater = distwater)

covariates_list <- lapply(covariates_list, function(im) {
  im$v[is.na(im$v)] <- median(as.vector(im$v), na.rm = TRUE)
  return(im)
})

model_null <- ppm(bee_ppp ~ 1, covariates = covariates_list)
print("=== NULL MODEL: ===")
## [1] "=== NULL MODEL: ==="
print(model_null)
## Stationary Poisson process
## Fitted to point pattern dataset 'bee_ppp'
## Intensity: 1.172671e-09
##              Estimate       S.E.   CI95.lo   CI95.hi Ztest      Zval
## log(lambda) -20.56398 0.02998801 -20.62276 -20.50521   *** -685.7402
model_proposed <- ppm(bee_ppp ~ 
                             elev + 
                             forest + I(forest^2) + 
                             hfi + I(hfi^2) + 
                             distwater,
                           covariates = covariates_list)
print("=== PROPOSED MODEL: ===")
## [1] "=== PROPOSED MODEL: ==="
print(model_proposed)
## Nonstationary Poisson process
## Fitted to point pattern dataset 'bee_ppp'
## 
## Log intensity:  ~elev + forest + I(forest^2) + hfi + I(hfi^2) + distwater
## 
## Fitted trend coefficients:
##   (Intercept)          elev        forest   I(forest^2)           hfi 
## -1.867799e+01 -3.617365e-03 -1.259259e-02 -1.137818e-05  4.956624e+00 
##      I(hfi^2)     distwater 
## -1.279527e+00 -1.427313e-04 
## 
##                  Estimate         S.E.       CI95.lo       CI95.hi Ztest
## (Intercept) -1.867799e+01 1.301992e-01 -1.893317e+01 -1.842280e+01   ***
## elev        -3.617365e-03 1.177883e-04 -3.848226e-03 -3.386504e-03   ***
## forest      -1.259259e-02 3.778325e-03 -1.999797e-02 -5.187206e-03   ***
## I(forest^2) -1.137818e-05 4.364851e-05 -9.692770e-05  7.417133e-05      
## hfi          4.956624e+00 5.414025e-01  3.895495e+00  6.017754e+00   ***
## I(hfi^2)    -1.279527e+00 5.869040e-01 -2.429837e+00 -1.292160e-01     *
## distwater   -1.427313e-04 2.046104e-05 -1.828342e-04 -1.026284e-04   ***
##                     Zval
## (Intercept) -143.4570263
## elev         -30.7107257
## forest        -3.3328493
## I(forest^2)   -0.2606775
## hfi            9.1551568
## I(hfi^2)      -2.1801295
## distwater     -6.9757578
print("=== AIC NULL MODEL: ===")
## [1] "=== AIC NULL MODEL: ==="
AIC(model_null)
## [1] 47960.3
print("=== AIC PROPOSED MODEL: ===")
## [1] "=== AIC PROPOSED MODEL: ==="
AIC(model_proposed)
## [1] 42064.23
print("=== ANOVA TEST: ===")
## [1] "=== ANOVA TEST: ==="
anova(model_null, model_proposed, test = "Chi")
## Analysis of Deviance Table
## 
## Model 1: ~1   Poisson
## Model 2: ~elev + forest + I(forest^2) + hfi + I(hfi^2) + distwater    Poisson
##   Npar Df Deviance  Pr(>Chi)    
## 1    1                          
## 2    7  6   5901.4 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Now we want to see whether the model is able to capture the correct shape:

intensity_im <- predict(model_proposed, type = "trend", n = 200)
log_intensity_im <- log(intensity_im)

plot(log_intensity_im,
     se = FALSE,
     superimpose = FALSE,
     main = "Log Fitted Intensity")

plot(log_intensity_im,
     se = FALSE,
     superimpose = FALSE,
     main = "Log Fitted Intensity")

plot(bee_ppp,
     pch = 16,
     cex = 0.6,
     cols = "green",
     add = TRUE)
## Warning in plot.ppp(bee_ppp, pch = 16, cex = 0.6, cols = "green", add = TRUE):
## 47 illegal points also plotted

med_elev <- median(elev, na.rm = TRUE)
med_forest <- median(forest, na.rm = TRUE)
med_hfi <- median(hfi, na.rm = TRUE)
med_dist <- median(distwater, na.rm = TRUE)


elev_effect <- effectfun(model_proposed, "elev", 
                         forest = med_forest,
                         hfi = med_hfi,
                         distwater = med_dist,
                         se.fit = TRUE)

forest_effect <- effectfun(model_proposed, "forest", 
                           elev = med_elev,
                           hfi = med_hfi,
                           distwater = med_dist,
                           se.fit = TRUE)

hfi_effect <- effectfun(model_proposed, "hfi", 
                        elev = med_elev,
                        forest = med_forest,
                        distwater = med_dist,
                        se.fit = TRUE)

dist_effect <- effectfun(model_proposed, "distwater", 
                         elev = med_elev,
                         forest = med_forest,
                         hfi = med_hfi,
                         se.fit = TRUE)

plot(elev_effect,
     legend = FALSE,
     main = "Effect of Elevation\nat median forest, hfi, distance")

plot(forest_effect,
     legend = FALSE,
     main = "Effect of Forest\nat median elev, hfi, distance")

plot(hfi_effect,
     legend = FALSE,
     main = "Effect of HFI\nat median elev, forest, distance")

plot(dist_effect,
     legend = FALSE,
     main = "Effect of Distance to Water\nat median elev, forest, hfi")

We further evaluate our proposed model using quadrat test, residuals, and partial residuals plots.

quadrat.test(model_proposed, nx = 5, ny = 5)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
## 
##  Chi-squared test of fitted Poisson model 'model_proposed' using quadrat
##  counts
## 
## data:  data from model_proposed
## X2 = 620.01, df = 14, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 21 tiles (irregular windows)
res <- residuals(model_proposed)
print(res)
## Scalar-valued measure
## Approximated by 3371 quadrature points
## window: polygonal boundary
## enclosing rectangle: [273874.9, 1870573.4] x [369042.8, 1735666.4] units
## 1112 atoms
## Total mass:
## discrete = 1112   continuous = -1112   total = -5.0343e-11
plot(res, main = "Model Residual", cols = "transparent")

par_res_elev <- parres(model_proposed, covariate = "elev")
par_res_forest <- parres(model_proposed, covariate = "forest")
par_res_hfi <- parres(model_proposed, covariate = "hfi")
par_res_dist <- parres(model_proposed, covariate = "distwater")

plot(par_res_elev,
     legend = FALSE,
     lwd = 2,
     main = "Partial Residuals: Elevation",
     xlab = "Elevation (m)")

plot(par_res_forest,
     legend = FALSE,
     lwd = 2,
     main = "Partial Residuals: Forest Cover",
     xlab = "Forest Cover (%)")

plot(par_res_hfi,
     legend = FALSE,
     lwd = 2,
     main = "Partial Residuals: HFI",
     xlab = "Human Footprint Index")

plot(par_res_dist,
     legend = FALSE,
     lwd = 2,
     main = "Partial Residuals: Distance to Water",
     xlab = "Distance to Water (m)")